Exact Solvability of the two-photon Rabi Hamiltonian 
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Exact spectrum of the two-photon Rabi Hamiltonian is found, proceeding in full analogy with the 
solution of standard (one-photon) Rabi Hamiltonian, published by Braak in Phys. Rev. Lett. 107, 
100401 (2011). The Hamiltonian is rewritten as a set of two differential equations. Symmetries that 
get hidden after further treatment are found. One can plainly see, how the Hilbert space splits into 
four disjunct subspaces, categorized by four values of the symmetry parameter c = ±1, ±i. There 
were only two values ±1 for the standard Rabi model. Four analytic functions are introduced by a 
recurrence scheme for the coefficients of their series expansion. All their roots yield the complete 
spectrum of the Hamiltonian. Eigenstates in Bargmann space are also at disposal. 

PACS numbers: 03.65.Ge, 02.30.1k, 42.50.Pq 



I. INTRODUCTION 

There were many trials to solve the model called (in 
quantum optics) Rabi Hamiltonian exactly, until Braak 
[l[ recently succeeded to do so. The same model is known 
under several pseudonyms, e. g. single-mode spin-boson 
system, Jaynes-Cummings model without rotating wave 
approximation and others; a brief survey can be found 
in 0]. Let us introduce a more general form of Rabi 
Hamiltonians Q, describing the interaction between a 
bosonic mode with energy u> and a two-level atom with 
level spacing ujq 



H 



(m) 



W 



+ utfb + g( ( T+ + a-) [(tf) m + b n 



(1) 



where m = 1,2,...; g is the interaction constant, a z , 
a ± are the Pauli matrices, and b are boson creation 
and annihilation operators, respectively. For the most 
studied m = 1 case Braak [l| recently presented an exact 
algebraic solution; we are going to solve the m = 2 model 
in full analogy. 

Before Braak's general solution, some special points, 
called also the Juddian ones, were known to be exactly 
solvable. This means that for some constraints on the 
Hamiltonian parameters, one can find both an eigen- 
function and its eigenenergy. It doesn't give the com- 
plete spectrum, just one excited eigenstate. This was 
shown at first for the standard m = 1 Rabi Hamiltonian 
Later the same picture, of course with different con- 
straints and eigenvalues was shown for m = 2, i. e. the 
two-photon Rabi Hamiltonian §. These special points 
manifest themselves as cross-sections of general solutions 
and they will serve as a check of our results. 

It is also known that the parameters of the two-photon 
Rabi model are restricted by \4g\ < uj, otherwise the 
eigenfunctions are not normalizable 
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oz 



->• z. 



(2) 



The eigenfunctions have to be analytic in the whole 
complex plane. Let us suppose that the solution of the 
stationary Schrddinger equation is a vector composed of 
two analytic functions (ipi(z),ip2(z)) T . We can insert 
eq.@ into the Hamiltonian (p}, case m = 2, and using 
the well-known 2x2 Pauli matrices, we obtain a coupled 
system of second-order differential equations 



2gip 2 ' + uzfa^ + 2gz 2 fa 
2gip'{ + ujztp' 2 + 2gz 2 fa - 



(«a _ E )fa = 



(*f + E)fa 



0, (3) 



where ' means the derivative with respect to z and E is 
the energy from the Hamiltonian's spectrum. For better 
symmctrization it is useful to go over to a linear com- 
bination of the considered functions, namely cf>i(z) = 
i>i (z) + ifaQz) and fciz) = ipiQz) - ifa(z)- Thus we get 
the set 



2g<p'{ + uzcj)^ 
-2g<j>' 2 ' + ujz<j)' 2 



(2gz* - 



E)<t>! + ^f0 2 = 
- E)<t> 2 + if fa = 0. 



(4) 



The next step is to find important symmetries present 
in the set (|4]) before they become less clear after some 
transformation followed by a series expansion. Here we 
prefer a somewhat different way than Braak [lj, but a 
very obvious one. We will perform two transformations 
of the variable z. The first one is z — > — z. One can easily 
see, that it leaves the set (j4j unchanged. Simple inspec- 
tion shows, that this implies common parity of both fa 
and fa. Plainly speaking, they must be either both even 



II. EXACT ALGEBRAIC SOLUTION 

The first step of the solution is going over to Bargmann 
space 0, introducing complex variable z, where the 
bosonic operators simplify to 



: 01 (2) 



(5) 



or alternatively, simultaneously odd 



2 



<j>i(-z) = -<j>i(z) 



(6) 



The second applied transformation is y — iz. The set 
(|3]) becomes 



-2fl 
2.9 ' 



dy 2 

2 02 



dy 2 



dy 2 

• + (2^-25)02 + ^*1=0. (7) 

dy 2 



Now the functions 0i and 02 evidently swapped their 
places, thus up to some common multiplicative constant: 



01 (iz) = c 2 (z) 

2 (iz) = c0i(z). (8) 

Possible values of the parameter 
inserting iz instead of z and we get <j 
c 2 4>i(z) and analogously (f> 2 (—z) = c 2 2 (z). Hence c 2 = 
1 from eq.([5| or c 2 = — 1 from eq.([n]). The symmetry 
parameter can acquire four values, c = ±1, ±i. In the 
standard m = 1 Rabi model such parameter had only two 
values ±1 and the transformation z — > — z was sufficient. 

Finally let us introduce four functions G c (z,E) whose 
roots with respect to E will be used to fulfill eq.®: 



c can be found by 
>i(~z) = c (j) 2 {iz) = 



G+(z) = 4> 2 (iz) - 0i(z) 
G-(z) = 2 (iz) + 0i(z) 
Gi(z) = i(f> 2 (iz) + 4>i(z) 
G-i(z) = i<f> 2 (iz) - 4>\{z). 



(9) 



In the last two cases the lower of eqs.© was multiplied 
by i in order to keep the G c functions real. They share 
the common parity of their 12 functions. The upper 
of eqs.((Hl) is then no more independent. The complete 
discrete spectra will be given by all roots of all G c (z, E) 
functions, again in full analogy with the standard Rabi 
model [lj]. Of course, the roots are meant with respect 
to E and they are independent of any chosen z. 

Having categorized the symmetries, we return to the 
set (|4| and perform the transformation 0i ( 2 = e~ KZ ipi,2- 
We get 



16ff 2 



8.9 



(11) 



In fact there should be ± in front of the square root, 
but the plus sign would make k divergent in the limit 
g 0, which is physically not reasonable. Notice that k 
remains real under the above mentioned restriction 4 \g\ < 
to. There is also a further analogy with the standard Rabi 
model, where the special case ojq = is exactly solved 
with help of the coherent state exp(±2g/u; 6')|0), here |0) 
is the lowest bosonic state 0, Q . In our notation, Braak 
performs the transformation 0^2 oc exp(— 2gz/w)0i i 2, re- 
call — > z in J2J). The m = 2 Rabi Hamiltonian solu- 
tion at uj = involves the squeezed vacuum term 
exp(±K 6^ 2 )|0). We will return to this case later. 

In the next step we expand the functions ^i,2 : 



n— — OG 
00 



(12) 



n— — 00 



To keep the solutions analytic we expect Q n {E) = 
K n (E) = for n < Q. Inserting these expressions 
into eas. (fT0|) we get the iteration scheme 



2g(n + 2)(n + l)Q n +2 + [(w - 85^)^ - Agn - E] Q r 



w 



K n = 



-2g(n + 2)(n + l)K n+2 + [(00 + 8gn)n + 4gn - E]K n 

-4uj K K n _ 2 + ^-Q n = 0. (13) 

Notice that the indexes differ by 0, 2 or 4, thus only the 
coefficients Q n , K n with common parity will be non-zero, 
either for even n = 0, 2, 4. . ., eq.©, or odd n = 1,3,5..., 
eq.©. The prefactor e~ KZ doesn't spoil the parity. The 
starting point of our iteration scheme is either at n = or 
n = 1 . Let us first have a look on the case n = and the 
symmetry parameter c = 1. Upper eq.® (f>i(iz) = 4> 2 {z) 
at z — implies Qq — K$, which can itself be a func- 
tion of Hamiltonian parameters, serving as normalization 
constant for the eigenfunctions. But as we are interested 
only in the roots of G+(z,E), this constant becomes an 
unimportant multiplication factor and we can choose 



2gi>'{ + (u>- 8gK)z{p[ - {4gn + £)<0i + <ftp 2 = 
-2g4)' 2 ' + (oj + 8gn)z-ip 2 - (Aujkz 2 - 4gn + E)^ 2 + 

+^1 = 0, (10) 

where we used the parameter k to simplify the first equa- 
tion by removing the term (8gn 2 — 2luk + 2c/)z 2 ?/5 1 , thus 
specifying its value 



Qo = 1, K Q = 1. 



(14) 



The second case is c = — 1. We have 0i(iz) = ~4> 2 (z) 
and again at z = 



Qo 



Kn 



(15) 
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The third case with c = i is only a bit more com- 
plicated. The symmetry condition from eq.@ is now 
4>i(iz) = i<p2{z) and its value at z = is rather trivial 
Qo = Kq = 0. But we will compare the first derivatives 
w. r. t. z, i. e. (f>i(iz) — i<j)' 2 {z) at z — and we get 
iQi = iKi, thus now our starting point is 



At last for 



Qi = 1, K l = 1. 

— i we get 



Qi = i, 



-1. 



(17) 



Concluding this part we can see that the Hilbert space 
of eigenf unctions splits into four disjunct subspaces. The 
corresponding eigenvalues can be found separately as 
roots of four G C {E) functions, eq.®. We substitute 

2 — — 

01,2 = e~ KZ ^1,2, the coefficients of expanded ipi,2 func- 
tions are found applying the iteration scheme (TIB")) subse- 
quently with four starting points, eqs. (|14IU"7"|) . The coef- 
ficients not defined by eq. (fl3|) are zero because of parity 
demands. K-2 = K-i = as well. 

Before proceeding to numerical calculations, we ought 
to mention several special cases, where the exact solution 
was already known. They will serve as a check of our 
general solution. 



III. SOME EXACTLY KNOWN CASES 

A. Case g = 

If the interaction constant is zero, the system sepa- 
rates into independent two-state atom with energy lev- 
els ±ujq/2 and a phonon with the modes Nuj, where 
N = 0, 1,2, . . . Thus the overall energy is ±uj /2 + Nuj. 
It is instructive to see how these values split into four 
groups as roots of four G c functions. Therefore we will 
solve this simple case explicitly. We return to the original 
functions "01,2, as the set of equations ([3]) decouples for 
.9 = and the two independent solutions are found easily 



'02 



2E4 

C 2 z— 



- = c lZ k ' 

- - r ? k 

— ^iz , 



(18) 



where we denoted the exponents by kl and k, as the par- 
ity conditions ([5]) and ([6]) are common for 01,2 and ipi,2 
functions. They force k' and k to be integer and the ana- 
lyticity demands make them non-negative. The energies 
E = — uj{)/2 + kuj and E = wo/2 + ft'w should be common, 
but they are generally different, thus the overall solutions 
will have either C\ = or Ci = 0. 

Let us first analyze the case (Q,ip2) T with the energy 
E = —ujq/2 + ku. The linear combinations 0i = "01 + 
•02 = C 2 z k and 02 = 4>i — -02 = —C2Z k . The solutions 
for the functions G c — finally yield 



G_ = -C 2 (iz) k + C 2 z k = = 
G+ = -C 2 {iz) k - C 2 z k = =j 
G- t - -iC 2 (iz) k + C 2 z k = =} 
d = -iC 2 {iz) k -C 2 z k = 



The second solution (?/>i,0) T with 0i = 02 = C\z 



• k = 0,4,8,... 

k = 2,6,10,... 

k = 3,7,11,... 
> k = 1,5,9,... (19) 

k' 



(16) and eigenenergies E = ojq/2 + k'ui gives 



Gi : 



Ci{iz) k +C lZ k =0 



= C 1 {izf -C x z k 
-- iCxiizf + C x z k 
iCxiizf -C x z k ' 



= 
= 
= 



k' = 2,6,10, 

• ft' = 0,4, 8, 

• k' = 1,5,9, 
k' = 3,7,11, 



(20) 



We will later analyze mainly the cases when loq and 
uj are comparable and the eigenenergies as roots of G c 
functions reorganize as follows: 



G_ : 


LU 


u> 


2 


2 


G+ : 






2 ' 


2 


G_, : 


Ll> 

2 


+ w, 


Gi : 


2 


UJ, - 



f 2u), 

wo 
2 

wo 
' 2 



— +4^, 



Wo 

2 

3w, 
3w, 



4ui, 



UJCl 

2 
2 



-~ + 5uj, 



5oj, 



(21) 



We can see that real values of the symmetry parameter 
c are connected with even number of phonon excitations 
N whereas the imaginary c is coupled to odd N. The 
global ground state is always given by the lowest root of 
G_, which is true also for non-zero g. 



B. 



Case loo = 



If the gap ujq between atomic levels disappears, all 
eigenenergies become twice degenerate. They are exactly 
known and now the complete spectrum is given by 



w , 1 , „ 



E = -- + (n + n = 0,1,2,... (22) 

where another dimensionless quantity was introduced 

8gn 



i 



16.g 2 



1 



(23) 



Let us reproduce this result. We return to the 0i.2 
functions, because now the system of equations (Q} de- 
couples for wq = and the general solutions are 



0i = exp 

Gii_ff_„_i (y^ff' 2 ) + G 2 i i-Fi {^p-, §, 
02 = exp (kz 2 ) • 



uifl J2 

ig z 



, (24) 
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exploiting the well known Hermitian polynomials H n and 
hypergeomctric function \F\. Further we introduced the 
quantity n = (bj+2E)/{2QiS) — l/2. It is nothing else but 
eq. (f2"2"j) reversed. Thus if the parity demands force n to be 
non-negative integer again, the spectrum is reproduced. 
Let us show it in detail at least for the simpler case of 
n even. We make use of a formula relating Hermitian 
polynomials and hypergeometric function, valid for n — 
0,2,4,... 



2»/2(n-l)!! 1 1 V 2 '2' 9 Z )' [ 

whereas for any n, including non- integer one, the Rum- 
mer transformation pdj gives 



1 1 



(26) 



First of all we set C\\ = 0, because the Hermitian 
polynomial with negative index has no parity, whereas 
the rest of the solutions (|2"4")l are even functions for n 



even. 
G± = 



We define q 2 — uifl/(4g), y = q 2 z 2 and require 



(27) 



using eqs. ([23 |) - (|26j) . We get 



{-2) n ' 2 {n - l)!!Ci 2 + G 22 T C21 = 0. 



For n non-integer (or odd integer) , the term with the Her- 
mitian polynomial would become complex, but our func- 
tions (j)2{iz) and G± are real - hence we set also Cyi = 0. 
We already mentioned that one parameter can be cho- 
sen, say the integration constant C21 = <^i(0) = 1. Thus 
C22 = <p2 (0) = ±1 and the two G c functions are iden- 
tical, G+ = G_. One can compare the series expansion 
of these (j^i^iz) solutions with those given by the scheme 
(TT3"1) . starting points (|T4")) or (|15l) and see that in fact we 
managed to perform the complete sum of eqs. (|12l) . 

Having exact formulas for some 01,2 and G c (z, E) func- 
tions at disposal allows us to make several observations. 
In the on-line supplement to ref. [1], Braak reports some 
problems with the radius of convergence R of his series 
in z. R seemed to be finite in some cases, though ana- 
lyticity in the whole complex plane of z is required. To 
keep the series analogous to our eqs. (|T2l convergent, as a 
necessary condition, the ratio K n+ i/K n had to go to zero 
for n — > 00. If it was non-zero, R became finite. For the 
special case ujq = we can calculate this ratio explicitly. 
The coefficients Q n with n = 2k become rather simple 



® 2k = (2 g )fc(2fc)' ( E ~ £2k ^ ( E ~ £2k ^ ■ ' ■ ^ E ~ ' 

(28) 

where e„ are eigenenergies from eq. (|22[) . En route we see 
that the energy E — e n terminates ipi (z) to a (Hermite) 
polynomial. The ratio 



^2fe+2 



IQ 



2k 



LOVL 1 

~2~5~2fc 



k —> 00 



(29) 



as required. The equations for K2k are not so simple, 
thus we resort to the exact solution of 4>2 and find 



K2k+2/K.. 



2k 



Wfi 1 

~2g2k, 



0, 



00. 



(30) 



For n odd and especially for non-zero wo w e performed 
at least numerical study of K n , Q n from the scheme fUS]) 
and found that their ratios are also proportional to \ jn 
in leading term. Hence we experience no problems with 
analyticity of the expanded functions including even very 
large z. 

The next remark concerns the practical numerical cal- 
culations of the roots G c (z,E) = 0, independent of z. 
It turns out that such calculations are numerically more 
stable for large z, which is allowed by the previous no- 
tion. One can exploit the large-z asymptotic for the <p± (z) 
contribution to G± given by [ll| 



1 F 1 (a,b,y) 



r(o) 



i-b 



Ol 1 - 

,V 



y >o, 

(31) 

which is analytic for non-negative integer power a — b = 
n/2, i. e. n = 0,2, ... as expected. <f>2(iz) gives the same 
after Kummer transformation (|26p . 

Further we return to the complete spectrum in eq. (|22p . 
For \g\ — > oj/A the quantity fl — > and the energy be- 
comes infinitely many times degenerate. This point is 
physically unsound, though well defined in the sense of a 
limit. 

Concluding, the eigenvalue of eq. (|22l) for n — coin- 
cides with mutually equal lowest roots of G_ and G+, 
for n — 1 it is the lowest root of both G_,; and Gi, for 
n = 2 the second lowest root of both G_ and G+, etc. 



C. Special cases with non-zero u)o,u,g 

Let us now recall the result of Emary and Bishop 0] , 
who found a set of isolated solutions for our model. We 
are not going to rederive it, but the basic fact is that 
under some constraint on Hamiltonian's parameters and 
energy E, at least one of the original eigenfunctions il>x,2 
becomes a product of some exponential function and of 
a polynomial. 

The main statement is that there exist exactly known 
eigenstates with eigenenergies 



E = — ^ + (n + ~ ) < L; 



(32) 



if the following conditions are fulfilled 



-0.5 1 1 1 1 1 1 1 1 1 *- 

0.05 0.1 0.15 0.2 0.25 

g 
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FIG. 1: Eigenenergies as roots of G C (E) functions for uio = 
and u = 1. Full lines are exact values. 



FIG. 2: Eigenenergies as roots of G C {E) functions for ojq — 1 
and u = 2. 



2 - 6S1 2 + —2- = 0, N = 2, 
4ur 



6 - ion 2 



4w 2 



0, N = 3, 



8(3 - 30ft 2 + 35ft 4 ) + 2(7 - 17ft 2 )-^% 



16cj 4 



0, A = 4, 



(33) 



etc. We will proceed so that we choose the values of 
ujq and u), than we gradually change g. Eq. ([33|) yields 
value(s) of g and followingly the energy is got from 
eg. (1321) . These solutions will manifest themselves as 
crossing points of appropriate roots of G_ and G+ in 
the case of N even, or crossings of and d roots 
in the case of A odd. This is again analogous with the 
standard m — 1 Rabi model, where the Juddian points 
appeared as crossings of the roots of G_ and 
Recall that contrary to eq. (|22p. this is only one known 
excited state for chosen set of Hamiltonian parameters, 
not the complete spectrum. 



IV. NUMERICAL RESULTS 

The resulting energies as roots of G c functions should 
not depend on z; this would be true after summing infi- 
nite number of summands in eq. (|12|) . In numerical cal- 
culations we truncate the sum and hope that the higher 
powers are not significant. We cannot use very small z, 
because of numerical instability. Somewhat surprisingly, 
we can use as large values as z = 1000 or even z = 10 4 
without real change of the roots. This fact was already 
noted for exactly solvable cases, nevertheless one even 
doesn't have to use compromise medium values despite 
of the truncated expansions. Using the symbolic program 
Mathematica, we could sum up to z L term with L = 34 



for even solutions or L = 35 for the odd ones. If the val- 
ues of a root with smaller L converged to the same value, 
we accept it. It turns out that for smaller values of \g\ 
the convergence is excellent, but it becomes poorer as we 
approach the maximal possible value, i. e. if \g\ — > uj/A. 
The interval of g in our plots is limited to < g < w / 4. 
It is known that although the eigenfunctions differ af- 
ter changing the sign of g, the eigenenergies remain the 
same, i. e. there is a mirror symmetry E{—g) = E{g). 

Let us test our calculations at first on the exactly 
solved case uo = 0. We choose u> = 1. Parts of parabolas 
with common top form the exact spectrum from eq. (|22l) , 
see full lines in Fig. 1. It is clear, that in the vicinity of 
the infinitely many times degenerate point with g = w/4, 
the roots of G C (E) functions become very dense and the 
functions themselves are quickly oscillating. That is the 
reason, why even for the orders as large as L — 34 the 
values of appropriate roots did not converge completely 
and we have to resort to some fitting procedure, yielding 
better guess of the saturation value for L — » oc. The 
roots are denoted so that the lowest one is Ei(G c ), the 
second lowest one E2(G C ), etc. We can see that the calcu- 
lated eigenenergies fit the exact values almost perfectly, 
except for some deviation at g close to w/4 and for the 
higher root, in this case E 2 (G±). 

In Fig. 2 we present three lowest roots of G±, full sym- 
bols, and two lowest roots of G±<, smaller empty sym- 
bols. The values u>o = 1 and ui = 2 are chosen so that 
the spectrum at g = becomes equidistant, see eq. (l2"Tl) . 
The large empty circles are exact solutions of eqs. (|3"3")) 
and (|32|) . The value of N is written nearby. We can see 
almost perfect match with the crossings of appropriate 
lines. The bottom four lines with N = and N = 1 do 
not cross; lines with N = 4 cross twice. 

Fig. 3 shows the same roots denoted by the same sym- 
bols as Fig. 2, just for ujq = 2 and uj = 1. Besides the well 
fitted exact crossing points, there are other crossings of 
lines with different N, which are not exactly known. A 
similar figure with in fact the same Hamiltonian param- 
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LU 




FIG. 3: Eigenenergies as roots of G C (E) functions for ujo — 2 
and uj = 1. The description of symbols is the same as in 
Figure 2. 



eters was already published @, where the authors plot 
also numerical results from larger matrices diagonaliza- 
tion. 

There is a couple of simpler analytic results, that can 
be got from our approach. So we can find e. g. the small- 
g expansion of the (global) ground state energy E = 
Ei(G-): 



+ 0{g 4 ), <7«w,w (34) 



2oj + ujq 

which can be compared with similar result for the m = 1 



Rabi model: « ~oj /2 - Ag 2 /(w + ujq) + . . . @. 



V. SUMMARY 

We have found the complete spectrum of the two- 
photon Rabi Hamiltonian as roots of four analytic func- 
tions G± and G±i, in the whole parametric space. These 
functions are given by the recurrence scheme (|13|) with 
four starting points (|14j) - (|17|) . The unnormalized eigen- 
functions in Bargmann space can be found as well, using 
%l>x(z) = [M*) + Mz)]/2 and ^ 2 (z) = - 2 (z)]/2. 

In his "Viewpoint: The dialogue between quantum 
light and matter" 12], E. Solano states that Braak [l[ 
managed to enlarge the class of exactly solvable models 
and that he added the standard Rabi model on the short 
list of exactly solvable quantum systems. We believe that 
this paper adds also the two-photon Rabi Hamiltonian on 
the same list. This list can almost surely be extended fur- 
ther by using Braak's approach on other related models 
and yet another task for future is deeper understanding 
of the criteria of its applicability. 
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